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We study the interaction driven localization transition, which a recent experiment in Gai-j^Mn^^As 
has shown to come along with multifractal behavior of the local density of states (LDoS) and the 
intriguing persistence of critical correlations close to the Fermi level. We show that the bulk of these 
phenomena can be understood within a Hartree-Fock treatment of disordered, Coulomb-interacting 
spinless fermions. The scaling analysis of the LDoS correlation demonstrates multifractality with 
correlation dimension d2 ~ 1.39. At the interaction-driven transition the states at the Fermi level 
become critical, whereas the bulk of the spectrum remains delocalized up to substantially stronger 
interactions. The mobility edge stays close to Fermi energy in a wide range of disorder strength, 
as the interaction strength is further increased. The localization transition is concomitant with the 
quantum-to-classical crossover in the shape of the pseudo-gap in the tunneling density of states, 
and with the proliferation of metastable HF solutions that suggest the onset of a glassy regime with 
poor screening properties. 



The recent tunneling experiments on Gai-^rMn^^As by 
Richardella et al. [.1^ have demonstrated critical, mul- 
tifractal correlations in the tunneling density of states 
at the metal-insulator transition in this semiconductor. 
While this would be expected from the theory of An- 
derson localization in non-interacting fermions, the ex- 
periment bears clear signs of the relevance of electron- 
electron interactions. Most strikingly, the critical corre- 
lations were found to persist very close to the Fermi level, 
even upon doping further into the metallic regime. This 
phenomenon clearly originates in interactions, which sin- 
gle out the Fermi level as a special energy. 

The effect of interaction has been studied on either 
side of the metal-insulator (MI) transition [2j , but rather 
little is known about its role close to criticality. In a 
seminal work Efros and Shklovskii (ES) showed [31 |3], 
that deep in the insulator the Coulomb interactions be- 
tween localized electrons create a pseudogap in the den- 
sity of states (DoS) near the Fermi level ep, where the 
DoS vanishes as p{e) ~ |e — e^P in 3D. This is reflected 
in the ES law of variable-range hopping conductivity at 
low temperatures, InfT(r) cx —^JTqJT. In the opposite 
limit of weakly disordered metals, Altshuler and Aronov 
(AA) discovered [5] interaction corrections to both the 
DoS near ep and the low T conductance. In particular, 
for spinless fermions, disorder and repulsive interactions 
both enhance the tendency to localize. In the weakly lo- 
calized regime the tunneling DoS has a dip at Ef with [S] 
(5p = p{£f + w) — p{£f) c<. \Ju}. However, in contrast 
to the classical ES pseudo-gap, the AA corrections are 
of purely quantum (exchange) origin: 8p cx f?!"^ , if the 
diffusion coefficient and the Fermi- velocity are held fixed. 

The quantum corrections of Ref. ^ can be effectively 
summed to obtain a non-perturbative result near d — 2 
by using the formalism of the non-linear sigma-model 
due to FinkeP stein [51 [7] . An effective action approach 



was suggested by Levitov and Sliytov ^ and Kamenev 
and Andreev to derive the non-perturbative expres- 
sion for the tunneling DoS in a weakly disordered 2D 
system near the Fermi energy. Remarkably, in the lower 
critical dimension d = 2, the DoS at sp vanishes ex- 
actly in the thermodynamic limit. In higher dimensions 
d > 2 instead, the above results suggest the following 
qualitative picture [UIZKII]: The pseudo-gap in the one- 
particle DoS gradually grows with increasing disorder or 
repulsion strength. p{sf) eventually vanishes at the lo- 
calization transition, and remains zero in the insulator. 
The shape of the pseudo-gap evolves from the quantum 
behavior 5p cx ^/uj in the metal to some non-trivial power 
p{ep -\- uj) (X ljJ^ , p > Q td, the Anderson transition point, 
to the classical behavior in the deep insulator. How- 
ever, the actual scenario might be more complex if the lo- 
calization transition is accompanied, or even preceded by 
a transition to glassy or other density-modulated phases. 

A fascinating property of electronic eigenfunctions 
near the localization transition of non-interacting parti- 
cles is their multifractality [12] . It is an exact property of 
critical states at the mobility edge e = but, also off- 
critical states exhibit multifractal character inside their 
localization or correlation radius ^ [13]. This fractality 
has important effects on interactions (both repulsive and 
attractive), as it enhances their local matrix elements. In 
the case of predominantly attractive interactions, it may 
induce local pairing gaps in weak Anderson insulators. 
In more conducting systems, it may lead to enhanced 
superconducting transition temperatures. |14H17| . 

It has remained an unresolved theoretical question 
whether such subtle wavefunction correlations survive in 
the presence of Coulomb interactions. On the other hand, 
the direct observation of multifractality in the tunneling 
spectra of Ref. [Tj strongly suggests so. Indeed, the mea- 
sured auto-correlation function of the local DoS (LDoS) 
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showed a well-established power-law decay with distance 
on the sample surface. Surprisingly, this critical behavior 
appeared nearly pinned to the Fermi energy without any 
fine-tuning of the Mn impurity concentration, implying 
that the mobility edge £m remains close to ei? in a broad 
range of disorder strengths. As mentioned above this in- 
dicates the importance of interactions, since they single 
out £f as the center of the pseudogap. In this Letter, 
we address the problem of multifractality at the inter- 
acting localization transition theoretically, and study the 
mechanism by which interactions pin the mobility edge 
Em nearly to ep in a broad parameter range. 

As a naive rationale for the pinning of the mobility 
edge one may consider that conditions for localization 
are more favorable where the density of states is low. 
Thus localization is naturally prone to occur first within 
the interaction-induced pseudogap, making £m track ep 
rather closely. However, it is only the single-particle (tun- 
neling) DoS that has a pronounced dip near ep, while 
the global thermodynamic DoS, dn/d^, that enters the 
conductivity via the Einstein relation, usually shows a 
different behavior. Hence, it is not obvious which notion 
of DoS is relevant for localization and transport purposes 
(see e.g. the discussion in [U [18] about global vs. local 
{dn/dij,)~^). A more detailed analysis is thus required in 
order to show that indeed the LDoS close to ep can be 
critical, while in the bulk of the spectrum the correlations 
are still metallic. 

To address these questions we consider a model of spin- 
less fermions on a 3D cubic lattice of size 10'^ with the 
tight-binding Hamiltonian 



(iJ) 



h.C. 



and interacting via long-range Coulomb repulsion: 

U sr-^ riirij 



Hi 



(1) 



(2) 



We employ periodic boundary conditions and choose 
t — 1 as the unit of energy. The on-site energies e.; 
are random, independently and uniformly distributed in 
ei £ [—^,^]- The chemical potential /i depends on in- 
teraction and is chosen so as to keep the average density 
1/2. For non-interacting particles {U — 0) the localiza- 
tion transition is known to occur at the disorder strength 
Wc = 16.5 in this model [H]. In the present work we 
choose W — U < Wc (not particularly close to Wc), 
so as to mimic conditions of Ref. ^ where the impurity 
concentration was not specially tuned. 

We attack this problem numerically by considering 
the interactions in the Hartree-Fock (HF) approxima- 
tion. This amounts to studying an effective single- 
particle model with self-consistent on-site energies and 
hopping amplitudes. We took into account the Hartree 
terms (occupation numbers nj in the sum rij) 



up to the 20*'' nearest neighbors, while the Fock terms 
(expectation values of cjcj) were considered only up to 
the fifth nearest neighbors [20]. Since the scaling prop- 
erties of critical states at the localization transition are 
independent of the details of the tight-binding Hamil- 
tonian, one might conclude that interactions in the HF 
approximation merely shift the critical point. However, 
this reasoning neglects the possibility of the hopping am- 
plitudes becoming long-ranged, as well as the emergence 
of long-range correlations in the disorder potential due 
to poorly screened Coulomb interactions. If those occur 
the existence of multifractality near the Anderson tran- 
sition is non-trivial, and the value of critical exponents 
may differ from the non-interacting case, even in the HF 
approximation. Our study indeed establishes the persis- 
tence of multifractality, albeit with slightly larger fractal 
dimension d2 ~ 1.39±0.05. This marks essential progress 
in comparison to earlier works [llj . As we will describe 
below, the critical behavior exhibits various further inter- 
esting features that are specific to the interacting case. 
Most importantly, we will establish that even in the in- 
sulating phase, the mobility edge remains very close to 
the Fermi level. 

-Pseudo-gap in the DoS. - One can easily show that 
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FIG. 1: (Color online) Disorder-averaged DoS p(e) of the HF 
states at W = 14 at different interaction strengths U. The 
crossover from the quantum AA correction 5p ~ -^/je — ef] 
to the classical ES gap p ~ (e — £f)'^ is seen. - Insert: p{£f) 
as a function of U. For U > 1.5 the DoS piep) ~ 0.5/{UL^) 
follows the classical ES law, where L = 10 is the system size. 

deep in the metallic and in the insulating regimes, HF 
correctly captures the AA and ES pseudogap features 
discussed above, while it provides a non-trivial mean 
field approach to describing various interesting phenom- 
ena happening at and close to the MI transition. In Fig.[l] 
we present the DoS of the HF levels. The curvature of 
the DoS p{e) at small uj = e — ep is seen to change sign as 
U increases. Theoretically one expects a critical power 
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law UJ^' for uj > l/p(w)C'^. We found [101 M « 0.6 ± 0-15 
in reasonable agreement with tunneling data in critically 
doped Si : B [U [3]. For smaller u, on the metallic side 
this turns into an AA dip Sp ~ y/uj, while an ES gap 
p ^ uj"^ is expected in the insulator [3 [TT] . 

In the insert of Fig. [l] we plot the DoS at the Fermi 
energy. For a finite sample size L well in the localized 
regime U > 1.5 it scales as p(e_F) ~ ijj2- A significant 
upturn sets in at L/ < 1, marking the crossover to quan- 
tum critical behavior near the localization transition |20j . 

-Auto-correlation of the local DoS. - To study multi- 
fractality of the local DoS, we have computed the spatial 
correlation of the HF wavefunctions |V'n(r)p: 



K{R;e) 



Er\Mr)\'\Mr + 'R)\'' 



(3) 



where e„ are the associated eigenvalues, (...) denotes the 
ensemble average over random realizations of on-site en- 
ergies e„, and is a narrow interval of energies of the 
order of the mean level spacing 5, centered at e. Mul- 
tifractal correlations [12l [15] imply that in the range of 
distances £o < R < the correlation function is the same 
as at the localization transition, K{R;e) ~ {£q/RY~'^^ . 
Here ^ is the localization or correlation length which di- 
verges at the transition, and £o is of the order of the 
lattice constant, d = 3 is the dimensionality of space and 
d2 < d is the correlation fractal dimension. For i? > ^ 
the correlation function K distinguishes delocalized and 
localized regimes, saturating to a constant in a metal and 
decreasing exponentially in an insulator. Close to criti- 
cality one expects scaling behavior, that is: the correla- 
tions should collapse to a single curve upon rescaling R 
by a (a finite size corrected version of ^), and amplitudes 
by /3, and expressing K{R;e) = fMj{R/ct), where 
fM,i{x) are universal scaling functions on the metallic 
and insulating sides of the transition, respectively. In or- 
der to optimize the choice of a, /3 (which both depend on 
U and e, while ly = 14 is fixed) and to determine the cor- 
relation dimension d2 wc use a simple analytical ansatz 
for /m,/, which captures the multifractal characteristics 
of the eigenfunction correlations: 



fi{x) = X 



(4) 
(5) 



The fits were optimized by the value i? sa 2. 

We first discuss K{R\£f) at the Fermi level, upon 
varying the strength of the interaction U , cf. Fig. [2j The 
good quality of the collapse demonstrates that the behav- 
ior of K{R; e) is consistent with multifractal correlations 
with dimension d2 ~ 1.39 ±0.05. This is only marginally 
larger than the one known for the non-interacting case 

d2 = 1.29 ±0.05 [ig. 

-MI transition. - Fig. [3] shows the evolution of 
the finite size-corrected correlation length a{U), as ob- 
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FIG. 2: (Color online) Data collapse of the auto-correlation 
of the LDoS, K{R;eF) onto (a) metallic and (b) insulating 
scaling functions /m,/ [Eqs. (|4|5[)]. We find d2 = 1.39 ± 0.05. 



tained from the scaling collapse. It exhibits strong non- 
monotonicity, indicating a localization transition: For 
U <U<^ 0.68, a{U) increases while for L/ > C/> « 0.74, 
it decreases. The best fits [20] to critical power laws yield 
^{U) = a\U ~U^^-^^\-•''''in with: 



U<i»\ 
a 0.50 ±0.05, 



« 1.08 ±0.05. 



(6) 



The difference in the fit exponents is too big to be a 
mere result of statistical errors, or of systematic errors 
related with the fitting procedure [50] ■ Indeed, as an 
independent check we computed the critical exponents 
for a non-interacting system of equal size, using the same 
method. We obtained [20l [22] W^^ = 19 ± 1, ^2 = 1-34 ± 
0.05 and much closer exponents I'm ~ 1-20 ± 0.05, lyj « 
1.08 ± 0.08. We thus believe that the difference in the 
fit exponents ([6]) is a genuine interaction effect. In this 
context it is interesting to note that the exponent um ~ 
0.5 has been reported in earlier experiments on Si : P, 
which remained a puzzle for theorists for a long time [7] . 

Eqs. (|6| might reflect the degradation of the multi- 
fractal pattern due to the interaction-induced mixing 
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FIG. 3: (Color online) Evolution of the finite size correlation 
length a with the interaction strength U at e — ef {W — U 
being fixed) . The raising and falling parts correspond to metal 
and insulator, resp. In the interval 0.7 < U < 0.8 the error- 
bars are too large to distinguish between various scenarios 
discussed in the main text. 

of non-interacting wavefunctions, which we expect to 
be much stronger in the delocalized than in the local- 
ized phase. Another phenomenon that undoubtedly in- 
fluences the MI transition is the gradual breakdown of 
screening in the metallic phase. The interactions, which 
are well screened deep in the metal, must become long 
range somewhere on the way to the insulator. This en- 
tails a crossover, or even a phase transition, to a glassy 
phase [5311^. We observe a trace of the latter via the on- 
set of non-uniqueness of the HF solutions roughly at the 
same point as the MI transition [25], but our resolution 
is not sufficient to determine whether the two phenom- 
ena coincide. It also remains an interesting open question 
whether screening breaks down at the MI transition only, 
or already within the metal, as it happens in mean field 
models with similar ingredients |26j . 

The large error bars in the interval U G [0.7, 0.8] do 
not allow us to determine ^(U) by an accurate treatment 
of the finite-size scaling a{U) = ^{U) fa{^/L). Two dif- 
ferent scenarii may be envisioned to reconcile the fit ex- 
ponents ^ with standard theoretical considerations: (a) 
there is a single localization transition close to Uc — 0.7 
with a shoulder on the metallic side due to an additional 
phase transition or a crossover in a different sector, such 
as the breakdown of screening or onset of glassiness. In 
that case vm would be expected to approach vi suffi- 
ciently close to Uc and on large scales, (b) there are 
two separate transitions at Ud ~ 0.68 and Uc2 ~ 0.74, 
with critical wavefunctions in an entire finite interval 
U G [Uci, Uc2]- In tbis more exotic scenario, there would 
be no a priori reason for the two exponents to coincide. 

Mobility edge in the bulk. - A similar scaling analysis 
as above may be performed for s away from sp 120' , which 



determines a critical line - the mobility edge £m{U) |27j . 
Our result is shown in Fig. |4] together with the band- 
edge, defined as the energy where p{e) drops to half of 
its maximal value. Fig. |4] demonstrates that the mobility 
edge is indeed trapped in a narrow range around sp, for 
U nearly all the way up to « 4 where the last states 
around the maximum of the DoS localize (at W = 14). 
This confirms the expectations of Ref. [1] that in a rela- 
tively broad region of parameters W and U, states near 
Ep are almost critical. Note also that is almost 5 times 
larger than U w 0.7 — 0.8 where the MI transition occurs. 
In fact U > U^, brings the system already very close to 
a Mott-type transition where charge-density wave order 
sets in. It remains an interesting question to study how 
such charge ordering effects and glassiness (i.e., the mul- 
tiplicity of HF solutions) affect the localization as the 
interaction strength increases. 
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FIG. 4: (Color online) Phase diagram. In a wide range of 
interaction strength \J the mobility edge (solid blue line) stays 
close to Ef as compared to the bandedges (green dashed line) . 

-Conclusion. - We have studied numerically the lo- 
calization transition in a 3D Anderson model of spinless 
fermions, with Coulomb interactions treated in HF ap- 
proximation. The MI transition was identified via the 
localization at the Fermi level, determined from a de- 
tailed study of the auto-correlation function of the HF 
eigenfunctions. Our main results are: (i) Multifractal 
power law scalings in the local DoS survive the presence 
of interactions, and extend up to a (large) correlation 
length £_{U,e). (ii) A significant critical Coulomb gap in 
the weakly insulating phase pins the mobility edge close 
to Ep for a wide range of parameters, while most of the 
higher energy excitations are still delocalized. In par- 
ticular, at disorder strength W — 14 (moderately close, 
but not fine-tuned, to the non-interacting critical disor- 
der W = 16.5) the critical Uc{ep) for the MI transition 
is ^ 5 times smaller than the C/* required for bulk local- 
ization. This is in good qualitative agreement with the 
experimental observations of Ref. [I] . (Hi) The apparent 
correlation length exponents display a significant asym- 
metry between the metallic and insulating sides, similar 
to tendencies reported in earlier experiments. We con- 
jecture that they arise from crossover phenomena in the 
metallic phase related with the breakdown of screening 
or the onset of glassy metastability scon in HF. Those 
deserve more thorough investigations in the future. 
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SUPLLEMENTARY MATERIAL 

I 



HARTREE-FOCK CALCULATIONS 

The effective Hartree-Fock (HF) Hamiltonian which 
corresponds to the model Eq. (1) in the main text, is: 

Hhf = ^ c\c, + i,, c\cj + h.c. . (7) 

i ij 

Here 



iij = -Uj - (c]c^)o, (9) 

where Uj = —t is the bare nearest-neighbor hopping 
and (...)o denotes the quantum-mechanical ground-state 
expectation value evaluated on the Slater determinant 
formed by the lowest N/2 HF levels. The effective on- 
site energy Vj contains the interaction-induced Hartree 
term which leads to correlated on-site energies, while the 
effective hopping tij contains the Fock term which may 
be long-range. In our calculations on a cubic 3D lat- 
tice of size L = 10 we take into account up to the 20**^ 
nearest neighbors (460 sites j nearest to i) in the Hartree 
term and the Fock terms corresponding to the 5^^ nearest 
neighbors (the 56 nearest sites up to distance \/5). We 
performed also a control calculation with both Hartree 
and Fock terms restricted to the 5*^ nearest neighbors, 
which confirmed all our main results close to the metal- 
insulator transition. Deep in the insulator side, it is how- 
ever important to keep the longer range Hartree terms to 
obtain the full classical Efros-Shklovskii gap, while long 
range Fock terms are negligible due to strong localization 
of the wavef unctions. 

The set of parameters X to be found self-consistently 
are all the (cjcj)o, (~ 28L^ parameters) plus the diag- 
onal parameters (rii)o = (c|^Ci)o. The chemical potential 
fi is always adjusted to assure half filling, as described 
below. 

In order to find the self-consistent solution we begin 
with a random initial guess for all the parameters X^^^ 
satisfying the condition 

^(n,)o=A/;, (10) 

i 

where the number of particles JVe ~ N/2 is fixed in 
our calculation (half- filling) . Diagonalizing the effective 

Hamiltonian Eqs. |7|-([9| using the initial x/"' one ob- 
tains the eigenfunctions ?/;,„ (r) and eigenvalues which 



allow us to compute the output parameters X^J^. : 

(n,)o = ^|V„(r,)P/(£m), (11) 
{4cj)o - 5]^:,(r.)V™(r,)/(£™)- (12) 

rn 

where /(e) is the Fermi distribution function with the 
chemical potential fi. It is to be found from the condition: 

M, = J2fiem). (13) 

m 

At r = considered in this paper there is an uncertainty 
of position of jj, between the two energy levels EA/'e+i ^^nd 
e jv^ . For most of the calculations we have chosen v — 
l(£A/'e+i + ^Af^)- However, for a more precise evaluation 
of the bottom of the DoS dip we used a parametric mixing 
/i = (1 — a) EA/'e+i + flSA^'c with < a < 1; a was fixed 
for a given disorder realization, but taken at random for 
different disorder realizations. 

An updated set of parameters to be used as initial pa- 
rameters X^^^^'^ for the next, i.e., (n + l)-th iteration is 
chosen as follows (n = 0, 1...): 

xl:+'^^il-a)X^:^+aX^:l. (14) 

The parameter a G [0, 1] is chosen such that the itera- 
tion process is stable and leads to a convergent solution. 
The iteration procedure is terminated and the output set 
of parameters is taken as the convergent solution if the 
absolute value of the difference between the values of all 
the parameters X of the previous and the final iteration 
is less than 10"**. 

Once the convergent solution is obtained and the final 
set of HF eigenfunctions '0m(r) and eigenvalues has 
been calculated, one can compute any quantity express- 
ible in terms of V'm(r) and The procedure is then 
repeated for different realizations of disorder to obtain 
disorder averaged quantities, such as the DoS and the 
LDoS correlation functions. 

DISORDER-AVERAGED DENSITY OF STATES 

The disorder averaged DoS is presented in Fig. [5] of 
the main text. Here we give a few more details of the re- 
sults. First of all we plot the Log-Log plot of the DoS at 
the Fermi energy, p{sf), as a function of the interaction 
strength. It is seen that in the insulating region, U > 1.5, 
piep) follows quite accurately the classical Coulomb glass 
law p{eF) = [7X2 with A « 0.51. However, closer to the 
localization transition at 0.75 < U < 1.5 p{ep) deviates 
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FIG. 5: DoS at the Fermi energy. For > 1.5 it is well 

A 



described by the "classical" expression p{£f) = ttt^ with 



A = 0.51 (the dashed line). In the delocalized phase U < 0.5 
the DoS is consistent with the critical behavior p{e; U) = 
p{e; 0) (1 - U/Ua)^ with 7 « 1.3. 



from the classical law, approaching quantum critical be- 
havior. In the delocalized phase it is consistent with 

p(ef)cx (1-C//C/,P, 7«1.3. (15) 

We also analyzed the behavior of p{e) close to £f- Close 
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FIG. 6: Log-Log plot of p(e) vs. {e — ef) in the metallic 
(17 = 0.25, U = 0.40 and U = 0.50), critical {U = 0.75), and 
insulating {U — 1.0, 1.5, 2.0) phases. This plot is given for 
comparison with the experimental results in Fig. 3 of Refs. [Ij 
and 131 

to criticality, where t = \1 — U/Uc\ ^ 1, one expects a 
scaling form for the density of states as 



with 



(16) 
(17) 
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FIG. 7: Collapse of data for the DoS in the window [7 3> a; S> 
5l (Eq. (27l) onto the scaling function fp{x): (a) metallic side 
with fp,M(x) = 1 -I- a;"'^^ (b) insulating side with fpj{x) = 
[x^^ -I- 



and 



fp{\x\ » 1) ^ \x\^ 



in the critical regime, whereas 
/p,M(|a;|<l) : 
/p,/(NI«l) ■ 



const., 

™2 



(18) 

(19) 
(20) 



which encompasses the above laws as limiting cases. 

The scaling is based on an assumption about the po- 
tential of a point charge within the critical regime, i.e., 
at a distance a <C r <C ^. Here ^ is the correlation length 
which diverges at the transition as. 



(21) 



and a is a certain microscopic length (e.g., the distance 
between donors in doped semiconductors [1]). In our sim- 
ulations it can be taken equal to the lattice spacing. The 
assumption is that the potential behaves as a modified 
power law 



V{r) - U 



(;) 



U - eVa. 



(22) 
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As Eq. ( 22 1 is essentially the relationship between the 



length scale r and the energy (or inverse-time) scale V, 
r] should coincide with the dynamical scaling exponent. 
For the non-interacting case the dynamic exponent takes 
its maximum value 77 = d = 3, while the minimal the- 
oretically admissible value cannot be smaller than the 
exponent of the Coulomb potential 0, ?7 > 1. 

The exponent 77 also governs the scaling of the static 
dielectric constant in the insulator 



(23) 



To show this it is enough to assume that in the insulator 
at distances r ^ ^ the potential V(r) takes the usual 
form of dielectric screening, 



Kor 



(24) 



and matches with the potential in Eq. ( 22 ) at distances 



The characteristic scale 6 in Eq. ( 16 ) is set by the po- 
tential V{r) at r ^ ^: 



(25) 



In a system of finite size L one should replace S by the 
mean level spacing Sl ^ V{L) in the critical region where 

Finally, the characteristic scale A of the DoS can be 
expressed through 5 and ^ by the obvious relationship: 



A" 



1 



U = e^/a. 



(26) 



Eqs. ( 25|26 l, as well as the scaling Eq. ( [T6| are valid for: 
a<<e, i7>w>(5L = C/(a/L)''. (27) 



From Eqs. ( |25|26| and ([Tt]) one immediately obtains the 
following scaling relations [1, 2 : 



7 = zy(3-r/), fi 



(28) 



in terms of rj and the correlation length exponent v. 

In Fig. [7] we verified the scahngs Eqs. ( 16|18|19l) b y 
collapsing the "low-energy" data (in the regime (|27[)) 
for p(e) close to onto the universal scaling functions 
fp.M{x) and fpj{x). From the power-law behavior of 
fp.M at X ^ 1 we found for the exponent /x: 



fiM = 0.53, 



(29) 



very close to the value experimentally observed in the 
tunneling DoS close to criticality |3]. We cross-checked 
this result in Fig. [S] by plotting In A vs. ln(l/(5). We ob- 
tained an almost linear curve in accordance with Eq. ( 17 1, 



with the same slope /i = 0.53 as in Fig. [7] In the insulator 




FIG. 8: Log- Log plot of A vs. 1/5, as obtained from the 
metallic side of the data collapse for the DoS. It is almost 
linear with the slope hm = 0.53. 



the best collapse corresponds to: 
m 0.68, 



(30) 



but the reliability of this exponent is less than the one 
in the metal (for instance a test similar to Fig. [s] yields a 
slope 0.54, consistent rather with ^m, but smaller than 
/i/ — 0.68). However, to be conservative we may conclude 
that the critical exponent /i lies in the range: 



/x = 0.60 ±0.15. 



(31) 



From the obtained exponent fi we can estimate the dy- 
namical scaling exponent rj: 



1 



= 1.9 ±0.2. 



A* 



(32) 



This gives the exponent ^ = 1.0 ± 0.3 in Eq. (23) in a 



reasonably good agreement with the experimental value 
l\ C « 0.71. 



LOCALIZATION TRANSITION AND 
MULTIFRACTAL CORRELATIONS 

Finite-size correlation length 

According to the phenomenology of multifractality the 
LDoS correlation function, normalized as in Eq. (3) of 
the main text, should have the form: 



(33) 



where f{x) is a scaling function and a and (3 are the 
finite-size correlation length and the finite-size correla- 
tion scale, respectively. For a; <C 1 the function f{x) is 
a non-trivial power-law a;"''^"''^-', where ^2 is the corre- 
lation fractal dimension, while at x ^ 1 it distinguishes 
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between the localized and delocalized regimes. In the de- 
localized regime f{x) = fnix) tends to a constant, while 
in the localized regime f{x) = fi{x) ^ e'^l^ decreases 
exponentially. In order to define a and /3 unambiguously 
we fix the asymptotic value of /m (2; — >■ 00) = 1 and the 
constant B in the exponent of //(x) to be the same as 
the one that controls the difference futix) — 1. 

The two scales a and /3 are related: /3 ~ {a/£o)^~'^^ , 
where £0 is the model-dependent ultra-violate cut-off of 
multifractality of the order of the elastic scattering mean 
free path. An important difference between the length 
scales a and £q is that the former is enhanced near the lo- 
calization transition, while the latter remains essentially 
constant. For an infinite sample the scale a coincides 
with the correlation/localization length ^, which diverges 
at the transition point. For a finite system of size L the 
divergence is rounded and the maximal value of a is of 
the order of L. Therefore we refer to the length scale 
a = a{U, W, e; L) as the finite-size correlation length. At 
a fixed system size this length depends on the parame- 
ters that control the localization transition: the disorder 
strength W, the strength of interaction U, and the energy 
e of the considered Hartree-Fock levels. 



Scaling collapse of data and the fractal dimension d2 

The goal of our study is two-fold: (i) to demonstrate 
that the above multifractal phenomenology does describe 
the data for the LDoS correlation function K{'R,e), in 
particular concerning the asymptotic power-law behavior 
of f{x), and to determine the correlation fractal dimen- 
sion ^2 ; and (ii) to identify the localization transition and 
find the corresponding critical behavior of the correlation 
length ^. To increase the accuracy and demonstrate scal- 
ing universality we collapse onto one single scaling curve 
/m or // all the data for K{Il, e) obtained by diagonal- 
ization of the HF Hamiltonian at different values of one 
of the control parameters t = W,U,e that drive the sys- 
tem across the localization transition. The quality of the 
collapse is monitored by the least deviation of all data 
points at a given value of r from the simple analytical 
interpolation formulae (Eqs. (4) or (5) of the main text) 
that capture the above multifractal phenomenology. In 
this way an optimal set of scaling parameters ^(t) and 
/3(t) is found for any value of t. Simultaneously, by fine- 
tuning the exponent 3 — ^2 of the power-law in Eqs. (4,5) 
and the number B in the exponential factor e^^/^ in 
fi{x) we found the values for the parameters c?2 and B 
that provide the best global fit. Note that a rough es- 
timation of the exponent ^2 can be obtained from the 
log-log plot of data points for one single value of r. The 
sensitivity of the quality of the collapse to the choice of B 
appeared to be very small, so that the rough estimation 
B = 2 was sufficient for our goals. 

The main result of the collapse is the precise determi- 



nation of (^2, obtained separately from the metallic and 
insulating sides of the transition. The resulting numbers 
showed relatively little fluctuations, even though they 
were slightly different in the metal and the insulator due 
to the finite range of t around the critical point. As 
a final estimate for d2 we take the arithmetic mean of 
the values obtained on either side of the transition. We 
emphasize that the resulting numbers do not take into 
account the finite-size corrections. This would require 
repeating the procedure for many system sizes, which is 
computationally very expensive. 

Using different parameters t to drive across criticality, 
we found the following values: 

d2 = 1.39(41, 37), £ = £f, W = U, t = U, (34) 
da = 1-35(26,44), C/ = 1.5, = 14, r = e, (35) 
d2 = 1-34(31,36), J7 = 3, 1^ = 14, r = e, (36) 

d2 = 1-34(31, 38), [/ = 0, e = 0, r = M/, (37) 

where the numbers in the parenthesis correspond to the 
two last digits in the metal and insulator fits. Note 
that restricting the Hartree terms from the 20"^ to 
the 5**^ nearest neighbors at e = ep changes ^2 only 
marginally ^2 = 1.38(38,37), approaching it towards the 
non-interacting case. The last set of data above corre- 
sponds to the non-interacting case treated by the same 
method for the same system size L = 10. This control 
measurement has been done for comparison, as the finite- 
size effects are presumably comparable in all these cases. 

We conclude that da stays pretty much constant inde- 
pendently of the choice of critical parameter and is very 
close to the non-interacting value. However, at the Fermi 
energy it is slightly larger than in other cases. Such an 
effect would not be entirely unexpected, since an analysis 
of localization in the presence of an interaction-induced 
pseudogap in lattices of high connectivity indeed sug- 
gested a suppression of the sparsity of wavefunctions at 
low energies. [1] 



Finite-size scaling for a{L) 

The finite-size correlation length a{L) can be generally 
represented as: 



a(L,r)=e(r)/„(e(T)/£), 



(38) 



where /^(a; — )■ 0) = 1 and fa{x^ 1) ~ 1/x. To fit data 
points corresponding to different values of the critical 
parameter r the correlation length ^ = ^(t) was chosen 
in the form: 



C = ii\l~r/T,\- 



(39) 
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where ii is of the order of the lattice spacing. 

Unfortunately, due to large error bars for a we were not 
able to determine the scaling function fa{x). Instead, we 
used the analytical fitting ansatz: 



a(L,r)=C(r) [I + C {ar) / Lf]-"^ , 



(40) 



which captures the main features of the supposed scaling 
function fa{x). It has the correct asymptotic limits at 
^ 3> i and ^ ^ L and has a regular expansion in t — r^, at 
i/^ ^ 1. The last condition applied separately for r < 
Tc< ("left" critical behavior) and r > rc> ("right" critical 
behavior) requires cr = l/v. However, if one assumes 
a single transition point at Tc with symmetric behavior 
of ^(r) around it, a regular expansion should start with 
(t — Tc)^. This requires a — Ijv. 

The parameters i\, Tc, C, and v were found by the best 
fit to the data air) separately on each side of the localiza- 
tion transition. It appears that except in the close vicin- 
ity of the transition where the error bars are very large, 
the finite-size corrections to a(V) are not very signifi- 
cant. Therefore, the apparent exponent v obtained from 



the power-law fit Eq. ( 40 1 in a region not very close to the 
transition point was largely insensitive to the finite-size 
corrections and could be obtained with reasonable accu- 



racy by neglecting the L-dependent term in Eq. ( 40 ) . 

Our fitting procedure consisted of three steps. In step 
1 we set C ~ and performed a three parameter fit for £i , 



Tc and v (see Eq. (39)) to roughly determine the critical 
exponent v. In step 2 we fixed the value of i> obtained 
in step 1, using it in Eqs. ( 39|40 | (with a ^ and 
a — 2/i/). We then made a three parameter fit for C, £o, 
Tc- Finally, in step 3, we fine-tuned the critical exponent 
1/ and the critical point Tc (as well as the length ii) by 
carying out a three parameter fit with fixed cr and C 
taken from the previous steps. The results are presented 
in Tables TIV, where the values characterize the last 
fine-timing step. 



TABLE I: Fit results for the interaction-driven localization 
transition at the Fermi energy e = ef, W — 14, for data 
obtained with Hartree terms up to the 20*'^ nearest neighbors, 
and Fock terms up to 5**^ nearest neighbors. 





a = l/u 


a = 2/v 


C = 




1.08 ±0.02 


1.05 ±0.03 


1.08 ±0.02 


I'M 


0.50 ±0.04 


0.47 ±0.04 


0.52 ±0.03 




0.68 


0.68 


0.70 


f/c> 


0.74 


0.73 


0.71 




0.23 


0.30 


0.23 


xli 


0.95 


0.93 


0.97 



TABLE II: Fit results for the interaction-driven localization 
transition at the Fermi energy e = ef, W = 14, for data ob- 
tained with Hartree and Fock terms up to 5**^ nearest neigh- 
bors. 





(7 = l/f 


a = 2/u 


C = 




1.31 ±0.11 


1.26 ±0.14 


1.31 ±0.11 




0.47 ±0.03 


0.43 ±0.03 


0.47 ±0.02 


t/c< 


0.76 


0.75 


0.78 


C/c> 


0.76 


0.75 


0.72 


X? 


0.69 


0.90 


0.69 


xlf 


1.14 


0.95 


1.09 



TABLE III: Fit results for the localization transition at the 
mobility edge e = £c away from the Fermi energy at fixed U = 
3, W = 14; with Hartree terms up to 20**^ nearest neighbors, 
Fock terms up to 5"^ nearest neighbors. 





a = l/u 


a = 2/u 


C = 




1.36 ± 0.12 


1.33 ± 0.12 


1.36 ±0.13 


I'M 


0.40 ±0.03 


0.38 ±0.13 


0.33 ±0.02 


£c< 


1.74 


1.75 


1.78 


£c> 


1.70 


1.70 


1.72 


X? 


0.19 


0.18 


0.19 


Xm 


0.10 


0.22 


0.11 



As a matter of fact, our data for a{T) do not allow the 
precise determination of the constant C, nor the crossover 



exponent cr. In view of the large uncertainty in C we used 
a fixed value of C = 0.5 in step 3, which is consistent with 
all the fitting attempts in step 2. However, the details of 
finite-size corrections do not change much the results for 
the correlation length exponent v on the metallic (I'm), 
and insulating (i^/) sides, nor do they affect much the 
critical point Tc obtained by the fit Eq. (39) from the 
sides of large t (rc>) and small r (tc<) sides. Within 
our accuracy they can be well determined even if the 
finite-size corrections ^/L are completely neglected (set- 
ting C = 0). Yet, the tests suggests that the fit with 
cr — Xjv is marginally better than the others. It is also 
the only reasonable fit from the theoretical viewpoint: 
On one hand, C = should be discarded on physical 
grounds. On the other hand, one would not expect the 
expansion of ^(t) to start with a (r — Tc)^ term (except 
at very small t — Tc which are beyond our reach) if the be- 
havior of a(T) is very asymmetric with respect to Tc. We 
will thus refer to the fits with a = l/v when discussing 
the general tendencies below. 

The first observation to make from Tables II-IV is that 
Vj grows with increasing interaction, while decreases 
as U grows. As a result, while the fit exponents are 
almost symmetric i^/ ~ ~1.15at?7 = 0, they become 
highly asymmetric i^j « 1.4, um ~ 0.4 at C/ = 3. 

A second observation from Table TII is that adding 
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TABLE IV: Fit results for the disorder-driven localization 
transition at e = for non-interacting particles. 





a = l/i^ 


a = 2/^ 


C = 




1.08 ±0.08 


1.06 ±0.07 


1.08 ±0.08 


I'M 


1.20 ±0.05 


1.16 ±0.05 


1.20 ±0.05 


Wc< 


19.8 


19.9 


20.3 


M^o 


17.9 


17.8 


17.7 


X? 


0.12 


0.13 


0.12 


Xm 


0.43 


0.48 


0.43 



long-range Hartree terms substantially increases uj while 
having a much smaller effect on vm- 

Both these observations are statistically significant. 

A third observation is statistically much less significant 
and concerns the sign of the difference Atc = Tc> — rc< 
between the apparent critical points obtained by fitting 
the data from the sides of larger and smaller values of the 
critical parameter r, respectively. Atc is found positive 
for the localization transition at the Fermi energy with 
long-range Hartree terms (Table I) but decreases sub- 
stantially when the Hartree terms are cut from the 20**^ 
down to the 5'^ nearest neighbors. Atc becomes nega- 
tive away from the Fermi energy (Table III) as well as for 
non- interacting particles (Table IV) . Note that this effect 
is not correlated with the asymmetry between and vm- 
Away from the Fermi energy (HI) the asymmetry is even 
larger than at £f (I), while in the non-interacting case 
(IV) it is much smaller than with interactions (I and HI). 
However, the difference Atc behaves similarly away from 
£f (HI) and in the non-interacting case (IV). Note that 
a statistically significant positive Atc would suggest the 
possibility of two transitions separated by a finite pa- 
rameter window with critical states. Such a possibility 
seems not to be excluded at the Fermi energy for suf- 
ficiently long-ranged interaction, but clearly requires a 
more refined numerical analysis. 

MULTIPLICITY OF HF SOLUTIONS AND 
GLASSY PHYSICS 

Finally we briefly address the issue of multiple solu- 
tions of the Hartree- Fock equations and its relation to the 
onset of glassy behavior as the interaction U increases. 

Our iterative procedure to solve the HF equations be- 
gins with an input configuration of occupation numbers 
nin(r) on each of the N lattice sites. At sufficiently large 
U the converged output HF solution nout(r) is generally 
different for runs with different inputs. To quantify this 
difference statistically we studied the quantity: 

^(t^) = ^ E E i"-t w - "iutWp, (41) 

m r 



where the superscript m denotes a solution obtained from 
a certain initialization n|™''(r), while denotes a refer- 
ence solution. In the simplest test we have chosen 10 
runs, of which 8 were with random inputs and 2 had a 
checkerboard order as input. The results for D{U) are 
presented in Fig.jOj One can see that for [/ < 4 the aver- 




FIG. 9: Average variance of on-site occupation numbers n(r) 
as a function of interaction strength. 
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FIG. 10: Difference between the output on-site occupation 
number and the checkerboard input occupation number: (a) 
17 = 4 output is random and uncorrelated with the input, (b) 
17 = 5 output has almost the same checkerboard structure as 
the input. 

age deviation of the solutions from the reference solution 
is small. It is thus reasonable to assume that physical 
properties evaluated on the various solutions are statis- 
tically very similar. However, starting from [/ « 4 the 
function D{U) sharply increases. 

In order to check whether this increase is due to a 
significant variation between random HF solutions or 
whether this increase in D{U) is due to stabilization of 
a checkerboard density pattern in the HF solution, we 
consider the solution obtained from initializing with a 
checkerboard input and plot the difference between the 
solution and the corresponding input pattern. 

The result for [/ = 4 shows (see Fig.[Io]^a)) that the dif- 
ference has a clear checkerboard structure which implies 
that the output was random. However, as U increases to 
U = 5 the difference reduces significantly (see Fig.[Io|^b)), 
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which signals the tendency to retain the checkerboard or- 
der in the solution. Comparing also with free energies of 
random solutions, we concluded that the transition to 
a checkerboard structure (charge density wave) occurs 
somewhere in the range 4 < [/ < 5. 

A more precise identification of the onset of multiplic- 
ity of solutions shows that it starts at much smaller values 
U « 0.7, which roughly coincides with the Uc at which 
localization at the Fermi energy occurs. In order to show 
this we generated 10 different HF solutions at t/ = 1.0 
and characterized them globally by the total energy E 
per site. The fact that the iterative HF procedure at 
U ~ 1.0 converges to different values of is a manifesta- 
tion of the existence of multiple local minima. Physically, 
one may expect that this will reflect in the onset of glassy 



behavior associated with the slow dynamics or relaxation 
between the minima that correspond to the various HF 
solutions. 

We then used the above solutions as inputs at slightly 
decreased interaction strength U — 0.9, the resulting so- 
lutions still being of different total energy (see Fig.jllja)). 
Upon decreasing the interaction in steps of 0.1, and using 
the solutions of the previous step as initial condition, we 
found that at J7 ~ 0.7 the total energies, after a large 
number of iterations, coincided (Fig. [Tl|(b)). That value 
of U can thus be interpreted, at this HF mean field level, 
as the border of a glassy regime. Upon further decrease 
of U the solutions did not diverge anymore, implying the 
existence of a unique HF solution (Fig. [TT|(c)). 




FIG. 11: Convergence of the HF procedure for the total en- 
ergy per site: (a) U — 0.9, the total energy converges to 
different output values, depending on the initial configura- 
tion of on-site occupation numbers; (b) U = 0.7, different 
initial configurations of occupation numbers obtained in the 
previous step converge to the same value of the total energy, 
and thus to the same HF solution; (c) U = 0.6 no further 
divergence in the total energy is observed any more. 
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